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Toward a universal water model: 

First principles simulations from the dimer to the liquid phase 
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A full-dimensional molecular model of water, HBB2-pol, derived entirely from "first principles" , is introduced 
and employed in computer simulations ranging from the dimer to the liquid. HBB2-pol provides excellent 
agreement with the measured second and third virial coefficients and, by construction, reproduces the dimer 
vibration-rotation tunneling spectrum. The model also predicts the relative energy differences between isomers 
of small water clusters within the accuracy of highly correlated electronic structure methods. Importantly, 
when combined with simulation methods that explicitly include zero-point energy and quantum thermal 
motion, HBB2-pol accurately describes both structural and dynamical properties of the liquid phase. The 
predictive power of the HBB2-pol quantum simulations opens the door to the long-sought molecular-level 
understanding of water under different conditions and in different environments. 



Given the central role that water plays in nature, it 
is not surprising that many studies have attempted to 
develop a microscopic picture of its unique properties. 
In particular, since the advent of computer simulations, 
myriad molecular models based on both force fields (in- 
cluding different degrees of empiricism) and ab initio ap- 
proaches have been proposed. Force field-based mod- 
els range from coarse-grained representations with no 
atomistic details (e.g, see Ref. 1), to classical param- 
eterizations in terms of point charges and rigid bonds 
(see Ref. 2) for a recent review), to yet more sophisti- 
cated models that account for molecular flexibility (e.g., 
Ref. |3|), electronic polarization (e.g., Ref. 0), and charge 
transfer- 5 . On the other hand, due to the associated com- 
putational cost, purely ab initio models so far have been 
limited to the use of density functional theory (DFT) 
(e.g., Ref. |i|). However, despite much recent progress, 
none of the existing models is capable of correctly de- 
scribing the properties of water from isolated molecules 
and small clusters in the gas phase up to the liquid and 
solid phases. Due to the broad spectrum of (often con- 
flicting) predictions derived from these models, the mi- 
croscopic behavior of water under different conditions 
and in different environments remains the subject of con- 
tinuing debate^"— . 

Independently of the number of molecules, the 
"first principles" modeling of water within the Born- 
Oppenheimer approximation builds upon two compo- 
nents: a faithful representation of the electronic poten- 
tial energy surface (PES) and the proper treatment of 
the nuclear motion at a quantum-mechanical level. In 
principle, the multidimensional PES can be accurately 
approximated at the coupled cluster level of theory in- 
cluding single, double and perturbative triple excitations, 
CCSD(T) 13 , which currently represents the gold stan- 
dard in quantum chemistry. For the nuclear degrees 
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FIG. 1. 2-body (a) and 3-body (b) interactions along with 
the second (c) and third (d) virial coefficients. Experimental 
data for second 33 and third 34 virial coefficients are shown as 
black squares and the HBB2-pol values as light blue circles. 



of freedom, quantum dynamics methods based on ba- 
sis set expansions of the nuclear wavefunction are well 
suited to the study of small complexes 14 , while simu- 
lation approaches based on the path-integral formalism 
allow the fully quantum-mechanical modeling of water 
in condensed phases^. In the end, these two compo- 
nents must be combined in an efficient computational 
scheme that enables the calculation of statistically con- 
verged quantities. 

Unfortunately, the computational cost associated with 
CCSD(T) makes these calculations prohibitively expen- 
sive even for small water clusters. To overcome this com- 
putational barrier while still providing an ab initio rep- 
resentation of the molecular interactions, less accurate 
but more efficient DFT approaches have been widely ap- 
plied to water simulations. However, the choice of the 
most appropriate DFT model for water remains the sub- 
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ject of ongoing researc h 16 ! 17 . Alternatively, the PES for 
a system of N water molecules can be expressed through 
the many-body expansion of interaction energies, consist- 
ing of 1-body, 2-body, . . . , N-body terms 18 . It has been 
shown that this expansion converges rapidly for water, 
such that it is sufficient to take into account only the first 
few 19 . Since the low-order interactions can be accurately 
calculated using CCSD(T), the many-body expansion ef- 
fectively enables the representation of the multidimen- 
sional PES at the CCSD(T) level of theory^. 

This strategy has been rigorously followed in the de- 
velopment of the CC-potfi and WHBB 29 models. In 
CC-pol, which describes the water molecules as rigid, 
the 2-body term was derived from CCSD(T) data while 
Hartree-Fock calculations were used to fit the 3-body 
term. WHBB, which allows for molecular flexibility, was 
parameterized using CCSD(T) and MP2 reference data 
for the 2-body and 3-body interaction terms, respectively. 
Both models include higher-body interactions through 
point polarizable dipoles. CC-pol accurately reproduce 
the vibration-rotation tunneling (VRT) spectrum of the 
water dimer and has been used in classical molecular dy- 
namics (MD) simulations of liquid water—. However, 
these simulations explicitly neglect nuclear quantum ef- 
fects such as zero-point energy and quantum thermal mo- 
tion, which have been shown to be important^. A new 
version of CC-pol with flexible monomers has recently 
been applied to study the water dimer 22 . WHBB also 
reproduces the VRT spectrum of the dimer with high 
accuracy 2 ^. Due to the computational cost associated 
with the high dimensionality of the polynomials used to 
represent the 3-body interactions, quantum simulations 
with the WHBB model so far have been limited to small 
clusters 24 . 

Here, we report on structural and dynamical proper- 
ties of water calculated with the newly developed full- 
dimensional model, HBB2-pol, derived entirely from first 
principles. Similarly to CC-pol and WHBB, HBB2-pol is 
built upon the many-body expansion of the molecular in- 
teractions. The 1-body term associated with intramolec- 
ular distortion is described by the spectroscopically accu- 
rate PES developed by Partridge and Schwenke 25 . The 
2-body interaction at short range is represented by the 
HBB2 potential 26 , which smoothly transitions between 
5.5 /AA and 7.5 /AA into the sum of electrostatic and 
dispersion interactions, reproducing the correct asymp- 
totic behavior. The induction contributions to non- 
pairwise additive interactions are taken into account us- 
ing Thole-type point polarizable dipoles on all atomic 
sites. In addition, an explicit 3-body component is in- 
troduced to account for short-range exchange-repulsion 
and charge transfer, which have been shown to make a 
significant contribution to the 3-body interactio n 27 ! 28 . 

The 3-body interaction of HBB2-pol is unique in that 
it is currently the first that has been fitted to CCSD(T) 
data. Importantly, the inclusion of induction interac- 
tions in the 3-body term enables the use of lower-degree 
polynomials than previously reported^ 9 -, resulting in a 
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FIG. 2. Relative energies of the low- lying isomers of the water 
tetramer (a), pentamer (b), and hexamer (c) computed with 
the HBB2-pol model (light blue circles) and state-of-the-art 
quantum chemistry methods (black squares) 35 i 36 . 



sizeable decrease in the computational cost associated 
with the many-body expansion. Specifically, the explicit 
3-body component takes the form of a third-degree poly- 
nomial in the exponential of the interatomic distances, 
which is invariant with respect to the permutations of 
equivalent atoms. This short-range correction was found 
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to be indispensable to overcome the limitations of current 
polarizable force fields (e.g., Ref. 1301 ), which attempt to 
describe all non-pairwise additive interactions through 
point polarizable dipoles. These findings suggest that 
the neglect of accurate short-range 3-body interactions 
invariably leads to an incorrect description of the liquid 
structure. 

The ability of HBB2-pol to accurately reproduce the 
molecular interactions of water is established by compari- 
son to high-level electronic structure calculations and ex- 
perimental measurements reported in the literature. Fig- 
ures and [T]3 show the correlation plots of the HBB2- 
pol and CCSD(T) 2-body and 3-body interaction ener- 
gies, respectively. For this analysis approximately 1400 
dimer and 500 trimer configurations were extracted from 
classical MD simulations of the hexamer, ice Ih, and liq- 
uid water. CCSD(T) energies were calculated with the 
aug-cc-pVTZ basis set 31 and corrected for the basis set 
superposition error through the counterpoise method 32 . 
Nearly perfect agreement is found over the entire range 
of energies, indicating that HBB2-pol accurately repro- 
duces the energetics of both favorably interacting and 
distorted dimers and trimers. 

An additional measure of the overall accuracy of the 
2-body and 3-body interactions is provided by the sec- 
ond and third virial coefficients, which directly probe 
dimer and trimer interactions, respectively. The HBB2- 
pol results are shown in Figures QJ and [Hi along with 
the corresponding experimental dat a 33 i 34 . In the calcu- 
lation of the second virial coefficient, nuclear quantum 
effects were explicitly included through the path-integral 
formalism (see Supporting Information for details). Since 
experimental data for the third virial coefficient are only 
available at relatively high temperatures where quantum 
effects are less important, these calculations were carried 
out at the classical level with each water monomer held 
fixed at the ground-state vibr at ionally- averaged geome- 
try. For both virial coefficients, the calculated values are 
in close agreement with the corresponding experimental 
data over the entire ranges of temperature, providing fur- 
ther support of the accuracy of the 2-body and 3-body 
interactions of HBB2-pol. 

The energetics of small water clusters with more than 
three molecules allows for a quantitative assessment of 
the ability of HBB2-pol to correctly describe N-body in- 
teractions with N > 3 and represents a stringent test 
for the overall many-body expansion. For this purpose, 
the relative energies of the low-lying isomers of the water 
tetramer, pentamer, and hexamer are compared in Fig- 
ure [2] with the corresponding ab initio values reported 
in the literatur e - 35 ! 36 . Water clusters in this size range 
are among the largest ones for which highly-correlated 
electronic structure calculations are still feasible. In all 
cases, HBB2-pol predicts the correct energy ordering of 
the isomers, with the energy differences being within the 
intrinsic accuracy of the ab initio methods 13 . 

The comparisons shown in Figures [T] and [2] demon- 
strate that the degree of accuracy of HBB2-pol is com- 
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FIG. 3. Oxygen-oxygen, 0-0 (a), oxygen- hydrogen, O-H (b), 
hydrogen- hydrogen, H-H (c), radial distribution functions. 
Experimental data from Refs. [37| and[38| are shown as dashed 
and solid black lines, respectively. The quantum RDFs calcu- 
lated with the HBB2-pol model are shown as light blue lines. 



parable to that of state-of-the-art quantum chemistry 
methods. While these highly-correlated calculations are 
restricted to small clusters, quantum simulations with 
HBB2-pol enable a microscopically detailed character- 
ization of water in condensed phases. To this end, 
HBB2-pol was used in fully quantum molecular dynam- 
ics simulations of liquid water at ambient conditions (T 
= 298.15 K and p = 0.997 g cm" 3 ). Specifically, path- 
integral molecular dynamics (PIMD) and centroid molec- 
ular dynamics (CMD) simulations were performed to de- 
termine structural and dynamical properties of the liq- 
uid phase, respectively^. Both PIMD and CMD are 
based on Feynman's formulation of statistical mechanics 
in terms of path-integrals and have been shown to ac- 
curately describe nuclear quantum effects in condensed 
phases of water—. All simulations were carried out with 
a system consisting of 256 molecules in a periodic cubic 
box. 

The oxygen-oxygen (O-O), oxygen- hydrogen (O-H), 
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TABLE I. Comparison between the experimental and cal- 
culated diffusion coefficient (D) and orientational relaxation 
time (Y2) of liquid water at ambient conditions. 





Experiment 


Simulation 


D (A'ps- 1 ) 


0.2^1 


0.23 it 0.5 


T 2 (ps) 


2.d!l 


2.5 ± 0.2 



a Ref. 40 
b Ref. El 



and hydrogen-hydrogen (H-H) radial distribution func- 
tions (RDFs) calculated from an 80 ps PIMD simulation 
are compared in Figure [3] with two sets of experimen- 
tal dat a 37 i 38 . The PIMD simulations correctly predict a 
lower first peak in the 0-0 RDF, as determined by the 
most recent experiments. This feature has been proven 
difficult to reproduce by current force field-based and ab 
initio models. Some differences between the PIMD re- 
sults and the experimental data exist for the second peak 
of the O-H RDF. However, both position and shape of 
this peak describing the spatial correlation between O 
and H atoms directly involved in hydrogen bonds are dif- 
ficult to be experimentally determined as demonstrated 
by the appreciable differences between the two sets of ex- 
perimental data. The diffusion coefficient (D) and orien- 
tational relaxation time (Y2) calculated at the quantum- 
mechanical level by averaging over six CMD trajectories 
of 10 ps each are listed in Table HI For both quantities, 
the CMD results are in quantitative agreement with the 
corresponding experimental values, providing further ev- 
idence of the accuracy of the HBB2-pol model. The aver- 
age structure of the first hydration shell was determined 
by labeling the water molecules according to the number 
of donating hydrogen bonds as non-donor, single-donor, 
and double-donor. Using the geometric criterion pro- 
posed in Ref-LD, the PIMD simulations predict that ~57% 
of the molecules are double donors and ~38% are single 
donors. Each molecule is involved, on average, in ~3 
hydrogen bonds. 

In summary, the full-dimensional HBB2-pol model, 
based entirely on first principles, has been introduced and 
employed in calculations of the water properties from the 
dimer to the liquid phase. Being derived from state-of- 
the-art quantum chemistry calculations, HBB2-pol rep- 
resents a major step toward the long-sought "universal 
model" capable of describing the behavior of water un- 
der different conditions and in different environments 39 . 
Future simulation studies with the HBB2-pol model will 
allow the resolution of current controversies regarding 
structural, thermodynamic, and dynamical properties of 
bulk, interfacial, and supercooled water. 
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